% load parameters (posterior mean) 


alpha1    =      1.79599541;  
alpha2    =      1.14322405;  
gamma1    =      0.09410810;  
gamma2    =      0.08092049;  
beta1       =      0.99914173;  
kappa1      =      0.44804350;  
tau1        =      3.98004910;  
rho_a      =      0.04901146;  
rho_f      =      0.94985429;  
rho_i      =      0.68761632;  
sd_a1   =      0.01528579;  
sd_f1   =      0.00325211;  
sd_i1   =      0.00614372;  
sd_a2   =      0.00693212;  
sd_f2   =      0.00115753;  
sd_i2   =      0.00324415;  
piess         =      0.00723860;  
gy         =      0.00496256; 
p11       =      0.96664397;  
p22       =      0.94965365;  
p11_f       =      0.90547855;  
p22_f       =      0.93909420;  

%--------------------------------------------------------------------
% default parameters
%--------------------------------------------------------------------


mhrund    = 1; 


%--------------------------------------------------------------------
% model descriptions and basic calculations
%--------------------------------------------------------------------

% call environments and model selection
rsmodel_spec2;

% parameter names
pnames = feval(fnames_.pnames,msel);


%--------------------------------------------------------------------
% means and standard deviations
%--------------------------------------------------------------------



P_i = [p11 1-p11;
       1-p22 p22];
   
P_f = [p11_f 1-p11_f;
       1-p22_f p22_f];   

P = kron(P_i,P_f);


%----------------------------------------------------------------------
% compute the solutions for inflation and output
%----------------------------------------------------------------------

  A_Aa_mat = [1-p11*rho_a, -(1-p11)*rho_a, -(p11*rho_a)/tau1, -((1-p11)*rho_a)/tau1, 1/tau1, 0;
        -(1-p22)*rho_a, 1-p22*rho_a, -((1-p22)*rho_a)/tau1, -(p22*rho_a)/tau1, 0, 1/tau1;
        -kappa1, 0, 1-beta1*p11*rho_a, -beta1*(1-p11)*rho_a, 0, 0;
        0, -kappa1, -beta1*(1-p22)*rho_a, 1-beta1*p22*rho_a, 0, 0;
        -gamma1, 0, -alpha1, 0, 1, 0;
        0, -gamma2, 0, -alpha2, 0, 1];

    solvec_a = inv(A_Aa_mat) * [rho_a/tau1;rho_a/tau1;0;0;0;0];

    A_Af_mat = [1-p11*rho_f, -(1-p11)*rho_f, -(p11*rho_f)/tau1, -((1-p11)*rho_f)/tau1, 1/tau1, 0;
        -(1-p22)*rho_f, 1-p22*rho_f, -((1-p22)*rho_f)/tau1, -(p22*rho_f)/tau1, 0, 1/tau1;
        -kappa1, 0, 1-beta1*p11*rho_f, -beta1*(1-p11)*rho_f, 0, 0;
        0, -kappa1, -beta1*(1-p22)*rho_f, 1-beta1*p22*rho_f, 0, 0;
        -gamma1, 0, -alpha1, 0, 1, 0;
        0, -gamma2, 0, -alpha2, 0, 1];
    solvec_f = inv(A_Af_mat) * [0;0;-kappa1;-kappa1;0;0];

    A_Ai_mat = [1-p11*rho_i, -(1-p11)*rho_i, -(p11*rho_i)/tau1, -((1-p11)*rho_i)/tau1, 1/tau1, 0;
        -(1-p22)*rho_i, 1-p22*rho_i, -((1-p22)*rho_i)/tau1, -(p22*rho_i)/tau1, 0, 1/tau1;
        -kappa1, 0, 1-beta1*p11*rho_i, -beta1*(1-p11)*rho_i, 0, 0;
        0, -kappa1, -beta1*(1-p22)*rho_i, 1-beta1*p22*rho_i, 0, 0;
        -gamma1, 0, -alpha1, 0, 1, 0;
        0, -gamma2, 0, -alpha2, 0, 1];
    solvec_i = inv(A_Ai_mat) * [0;0;0;0;1;1];

   
    pi1a = solvec_a(3);  pi2a=solvec_a(4);
    pi1f = solvec_f(3);  pi2f=solvec_f(4);
    pi1i = solvec_i(3);  pi2i=solvec_i(4); 
    
    i1a  = solvec_a(5);  i2a = solvec_a(6);
    i1f  = solvec_f(5);  i2f = solvec_f(6);
    i1i  = solvec_i(5);  i2i = solvec_i(6); 
   
Vs1 = diag([sd_a1^2/(1-rho_a^2);sd_f1^2/(1-rho_f^2);sd_i1^2/(1-rho_i^2)]);
Vs2 = diag([sd_a2^2/(1-rho_a^2);sd_f2^2/(1-rho_f^2);sd_i2^2/(1-rho_i^2)]);
    
   
   disp('relative volatility of interest rate: active, high vol regime'); 
   sqrt([i1a i1f i1i]*Vs1*[i1a i1f i1i]')/sqrt([pi1a pi1f pi1i]*Vs1*[pi1a pi1f pi1i]')
    
   disp('relative volatility of interest rate: active, low vol regime'); 
   sqrt([i1a i1f i1i]*Vs2*[i1a i1f i1i]')/sqrt([pi1a pi1f pi1i]*Vs2*[pi1a pi1f pi1i]')

  disp('relative volatility of interest rate: passive, high vol regime'); 
   sqrt([i2a i2f i2i]*Vs1*[i2a i2f i2i]')/sqrt([pi2a pi2f pi2i]*Vs1*[pi2a pi2f pi2i]')
    
   disp('relative volatility of interest rate: passive, low vol regime'); 
   sqrt([i2a i2f i2i]*Vs2*[i2a i2f i2i]')/sqrt([pi2a pi2f pi2i]*Vs2*[pi2a pi2f pi2i]')
			